clear all;
% clc;
warning off;

patnum=2;
ptag=['p',num2str(patnum),'_'];

diam   = [3,6,9,12,15];
diamsw = [3,6,9,12,15];
nodiam = length(diam);
nodiamsw = length(diamsw);

noaxons = 200;

epi_dcdata = zeros(noaxons,nodiam);
epi_drdata = zeros(noaxons,nodiam);
sub_dcdata = zeros(noaxons,nodiam);
sub_drdata = zeros(noaxons,nodiam);

epi_dcdatasw = zeros(noaxons,nodiamsw);
epi_drdatasw = zeros(noaxons,nodiamsw);
sub_dcdatasw = zeros(noaxons,nodiamsw);
sub_drdatasw = zeros(noaxons,nodiamsw);

% ----------------------- Additional Data ---------------------------------

iunit_epi = [1.8,1.5,1.6,1.8,1.7]*1e-3;
iunit_sub = [1.3,1.4,1.2,1.1,1.3]*1e-2;

ra_epi = 2./iunit_epi/1e3;
ra_sub = 2./iunit_sub/1e3;

isens_expsub = [0.20,0.35,0.70,0.55,0.80];
ipain_expsub = [0.50,0.50,1.50,1.20,3.10];

isens_expepi = [1.30,19.0,7.00,3.00,9.00];
ipain_expepi = [2.20,24.5,19.0,5.00,19.0];    

xo = 0;
xf = 12;
x = xo:(xf-xo)/(nodiam-1):xf;
isens_epimin = min(isens_expepi)*ones(size(x));
isens_epimax = max(isens_expepi)*ones(size(x));
isens_submin = min(isens_expsub)*ones(size(x));
isens_submax = max(isens_expsub)*ones(size(x));

ipain_epimin = min(ipain_expepi)*ones(size(x));
ipain_epimax = max(ipain_expepi)*ones(size(x));
ipain_submin = min(ipain_expsub)*ones(size(x));
ipain_submax = max(ipain_expsub)*ones(size(x));

isens_epiavg = median(isens_expepi)*ones(size(x));
isens_subavg = median(isens_expsub)*ones(size(x));
ipain_epiavg = median(ipain_expepi)*ones(size(x));
ipain_subavg = median(ipain_expsub)*ones(size(x));


% ----------------------- Import Data -------------------------------------

box_groups   = {'3','6','9','12','15'};
box_groupssw = {'6','9','12','15'};

cd([ptag,'mrgpotentials']);
model_iunit_epi=load([ptag,'epi_bp_current.txt']);
model_iunit_sub=load([ptag,'sub_bp_current.txt']);
cd ..;

cd([ptag,'mrgthresholds\']);
for ii = 1:nodiam
    
    dtag = [num2str(diam(ii)),'um'];
    temp_a = load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
    temp_b = load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
    temp_c = load([ptag,'sub_bp_dc_act_',dtag,'.txt']);
    temp_d = load([ptag,'sub_bp_dr_act_',dtag,'.txt']);
    epi_dcdata(:,ii) = abs( temp_a(:,1) )/(2/model_iunit_epi/1e3);
    epi_drdata(:,ii) = abs( temp_b(:,1) )/(2/model_iunit_epi/1e3);
    sub_dcdata(:,ii) = abs( temp_c(:,1) )/(2/model_iunit_sub/1e3);
    sub_drdata(:,ii) = abs( temp_d(:,1) )/(2/model_iunit_sub/1e3);
    
end
cd ..;

x_diam = ( diam'*ones(1,noaxons) )';

cd([ptag,'swthresholds\']);
for ii = 1:nodiamsw;
    
    dtag = [num2str(diamsw(ii)),'um'];
    temp_a = load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
    temp_b = load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
    temp_c = load([ptag,'sub_bp_dc_act_',dtag,'.txt']);
    temp_d = load([ptag,'sub_bp_dc_act_',dtag,'.txt']);
    epi_dcdatasw(:,ii) = abs( temp_a(:,1) )/(2/model_iunit_epi/1e3);
    epi_drdatasw(:,ii) = abs( temp_b(:,1) )/(2/model_iunit_epi/1e3);
    sub_dcdatasw(:,ii) = abs( temp_c(:,1) )/(2/model_iunit_sub/1e3);
    sub_drdatasw(:,ii) = abs( temp_d(:,1) )/(2/model_iunit_sub/1e3);
    
end
cd ..;

x_diamsw = ( diamsw'*ones(1,noaxons) )';

% ----------------------- Plotting DC Fibers ------------------------------

boxwid=1;
boxspc1=0.05; % spacing between a pair
boxspc2=0.15; % spacing between pairs
boxip=1; % initial position
posbox=zeros(1,10);
posbox(1)=boxip;
for ii=2:10
    if(mod(ii,2)==0)% even
        posbox(ii)=boxip+(ii-1)*boxwid+(ii/2)*boxspc1+(ii/2-1)*boxspc2;
    else% odd
        posbox(ii)=boxip+(ii-1)*boxwid+((ii-1)/2)*(boxspc1+boxspc2);
    end
end

xo_epi = posbox(1)-boxwid*(3/4);
yo_epi = isens_expepi(patnum);%isens_epimin(1);
xf_epi = posbox(end)+boxwid*(3/4);
yf_epi = ipain_expepi(patnum);%isens_epimax(1);
ym_epi = median(isens_expepi);

xo_sub = posbox(1)-boxwid*(3/4);
yo_sub = isens_expsub(patnum);%isens_submin(1);
xf_sub = posbox(end)+boxwid*(3/4);
yf_sub = ipain_expsub(patnum);%isens_submax(1);
ym_sub = median(isens_expsub);

atrans = 0.9;
xtloc = (posbox(1:2:9)+posbox(2:2:10))/2;

% epidural
figure;
hold on;
hr_epi=patch([xo_epi,xf_epi,xf_epi,xo_epi],...
             [yo_epi,yo_epi,yf_epi,yf_epi],'g');
alpha(hr_epi,atrans);
data = [epi_dcdata,epi_dcdatasw];
diamlab = repmat({'3' '6' '9' '12' '15'},1,2);
simobs = [repmat({'MRG'},1,5),repmat({'SW'},1,5)];
h = boxplot(data,{diamlab,simobs},'colors',repmat('kr',1,5),...
    'whisker',100,'factorseparator',1,'labelverbosity','minor',...
    'positions',posbox,'widths',boxwid*ones(1,10));    
hfp = findobj(gca,'Color',[0.75,0.75,0.75]);
hmrg = findobj(gca,'Color','k');
hswe = findobj(gca,'Color','r');
set(hmrg,'LineWidth',4,'LineStyle','-');
set(hswe,'LineWidth',4,'LineStyle','--','Color','k');
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',xtloc,'xticklabel',box_groups,'FontSize',36); 
set(hfp,'Color','k','LineWidth',2,'LineStyle','-');
hold off;
% title('Epidural DC','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',36);
ylabel({'DC Stimulation';'Threshold Current (mA)'},'FontSize',36);
xlim([xo_epi,xf_epi]);
ylim([0.5,1500]);

% intradural
figure;
hold on;
hr_sub=patch([xo_sub,xf_sub,xf_sub,xo_sub],...
             [yo_sub,yo_sub,yf_sub,yf_sub],'g');
alpha(hr_sub,atrans);
data = [sub_dcdata,sub_dcdatasw];
diamlab = repmat({'3' '6' '9' '12' '15'},1,2);
simobs = [repmat({'MRG'},1,5),repmat({'SW'},1,5)];
h = boxplot(data,{diamlab,simobs},'colors',repmat('kr',1,5),...
    'whisker',100,'factorseparator',1,'labelverbosity','minor',...
    'positions',posbox,'widths',boxwid*ones(1,10));
hfp = findobj(gca,'Color',[0.75,0.75,0.75]);
hmrg = findobj(gca,'Color','k');
hswe = findobj(gca,'Color','r');
set(hmrg,'LineWidth',4,'LineStyle','-');
set(hswe,'LineWidth',4,'LineStyle','--','Color','k');
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',xtloc,'xticklabel',box_groups,'FontSize',36); 
set(hfp,'Color','k','LineWidth',2,'LineStyle','-');
hold off;
% title('Intradural DC','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',36);
ylabel({'DC Stimulation';'Threshold Current (mA)'},'FontSize',36);
xlim([xo_sub,xf_sub]);
ylim([0.09,300]);

return

figure;
subplot(1,2,1);
hold on;
plot(x,isens_epiavg,'k--','LineWidth',4);
plot(x,ipain_epiavg,'k-','LineWidth',4);

hb1 = boxplot(epi_dcdata,'labels',{'','','','',''},'whisker',100,...
              'plotstyle','compact','color','k');   
hb2 = boxplot(epi_dcdatasw,'labels',box_groups,'whisker',100,...
             'symbol','k+','outliersize',25,'boxstyle','outline','colors','b');   
         
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);       
hold off;
title('Epidural DC','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',30);

set(hb1(1,:),'LineWidth',5);
set(hb1(2,:),'LineWidth',30);
set(hb1(3,:),'Marker','s','MarkerSize',10);
set(hb1(4,:),'MarkerSize',1);
set(hb2,'LineWidth',3);
legend('Sensory Threshold','Pain Threshold','Locatin','NE');
ylim([0.5,300]);

subplot(1,2,2);
hold on;
plot(x,isens_subavg,'k--','LineWidth',4);
plot(x,ipain_subavg,'k-','LineWidth',4);

hb1 = boxplot(sub_dcdata,'labels',{'','','','',''},'whisker',100,...
              'plotstyle','compact','color','k');   
hb2 = boxplot(sub_dcdatasw,'labels',box_groups,'whisker',100,...
             'symbol','k+','outliersize',25,'boxstyle','outline','colors','b');   
         
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);       
hold off;
title('Intradural DC','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',30);
set(hb1(1,:),'LineWidth',5);
set(hb1(2,:),'LineWidth',30);
set(hb1(3,:),'Marker','s','MarkerSize',10);
set(hb1(4,:),'MarkerSize',1);
set(hb2,'LineWidth',3);
ylim([0.1,50]);

return

% ----------------------- Plotting Intradural Case ------------------------

figure;
subplot(1,2,1);
hold on;
plot(x,isens_epiavg,'k--','LineWidth',4);
plot(x,ipain_epiavg,'k-','LineWidth',4);

hb1 = boxplot(epi_drdata,'labels',{'','','','',''},'whisker',100,...
              'plotstyle','compact','color','k');   
hb2 = boxplot(epi_drdatasw,'labels',box_groups,'whisker',100,...
             'symbol','k+','outliersize',25,'boxstyle','outline','colors','b');   
         
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);       
hold off;
title('Epidural DR','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',30);
set(hb1,'LineWidth',3);
set(hb2,'LineWidth',3);
ylim([0.5,300]);

subplot(1,2,2);
hold on;
plot(x,isens_subavg,'k--','LineWidth',4);
plot(x,ipain_subavg,'k-','LineWidth',4);

hb1 = boxplot(sub_drdata,'labels',{'','','','',''},'whisker',100,...
             'symbol','k.','outliersize',25,'boxstyle','outline','colors','k','notch','on');
   
       
hb2 = boxplot(sub_drdatasw,'labels',box_groups,'whisker',100,...
             'symbol','k+','outliersize',25,'boxstyle','outline','colors','b');   
         
set(gca,'FontSize',26,'YScale','Log','LineWidth',2);
set(gca,'xtick',1:size(box_groups,2),'xticklabel',box_groups,'FontSize',26);       
hold off;
title('Intradural DR','FontSize',30);
xlabel({'';'Fiber Diameter (\mum)'},'FontSize',30);
ylabel('Stimulation Threshold Current (mA)','FontSize',30);
set(hb1,'LineWidth',3);
set(hb2,'LineWidth',3);
ylim([0.1,50]);

